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Abstract 


Adjoint methods suitable for obtaining sensitivity derivatives for numerical simulations of solid-oxide fuel cells are presented. The adjoint 
method is derived, and the implementation is discussed, including a methodology for accurately obtaining all the linearizations necessary for 
correct implementation. Results are included for a one-dimensional anode model that includes diffusion, permeation, and relevant chemical 
reactions. Using this model, the accuracy of the sensitivity derivatives is demonstrated for design variables describing geometric and material 
properties of the anode. Finally, the adjoint method is demonstrated for a three-dimensional fuel cell geometry where sensitivity derivatives are 
obtained for approximately 180,000 design variables. The results are used to modify the upper and lower walls of the plenum to obtain significantly 


improved distribution of fluid amongst the channels. 
© 2007 Elsevier B.V. All rights reserved. 
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1. Introduction 


Numerical simulations of solid-oxide fuel cells (SOFC) have 
been used to gain understanding of important physical phenom- 
ena and to supply guidance in the continuing development of 
improved fuel cells [1-8]. To date, the simulations have been 
primarily focused on analysis of fuel cells or fuel cell compo- 
nents, without strong emphasis on utilizing the simulations in a 
design optimization environment. Because of the emphasis on 
analysis instead of design, sensitivity information to determine 
the effects of variations in design parameters on performance 
has been primarily implemented by simply changing the param- 
eter of interest, re-running the simulation, and comparing the 
results with those from the original simulation [1,5,6,8]. While 
this approach can be used to determine the effects of parameter 
variations on fuel cell performance, a more rigorous approach 
toward optimization would likely lead to better designs, and 
can also provide improved insight into the parameters affecting 
the performance of the fuel cell. For SOFC problems, exam- 
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ple cost functions that can be used for improving performance 
include minimizing temperature variations, obtaining equal dis- 
tribution of fuel in each of the channels, or maximizing power. 
Although not included in this study, time-dependent formula- 
tions are also possible that may be useful for minimizing start-up 
and short-term transient times. Design variables may be related 
to the shape/size of the fuel channels, electrodes, electrolyte, 
and interconnect, but may also be coupled to the stoichiometric 
composition of fuel or material properties such as the porosity 
or tortuosity of the electrodes. 

In refs. [9,10], optimization algorithms have been used to 
improve the performance of a polymer-electrolyte-membrane 
fuel cell (PEM) using four design variables, where the sen- 
sitivity derivatives used for the optimization algorithm have 
been obtained using a finite-difference approach. While finite 
differences are often a viable means for computing sensitiv- 
ity derivatives, this method can be computationally restrictive 
when a sufficiently large number of design variables are present. 
In addition, accurate derivatives can sometimes be difficult to 
obtain using finite differences because of subtractive cancella- 
tion errors [11,12], which occur when the function evaluations in 
the numerator become computationally indistinguishable when 
very small perturbations are used. To date, there has not been 
extensive research targeted at providing accurate sensitivity 
derivatives that may be used in conjunction with optimization 
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Nomenclature 

B effective permeability coefficient (m? s7!) 

c molar concentration vector (kmol m~?) 

Dij binary diffusion coefficient (m? s7!) 

D} effective diffusion coefficient (m? s7!) 

Dp Knudsen diffusion coefficient (m? s7!) 

f cost function (cost function dependent) 

h mesh size (m) 

I current density (A m~?) 

J mass flux vector (kg m~? s7!) 

Kn Knudsen number 

kb backward reaction rate (reaction type dependent) 

kf forward reaction rate (reaction type dependent) 

L augmented cost function (cost function depen- 
dent) 

M molecular weight (kg kmol~!) 

N molar flux vector (kmol m~? s7!) 

p pressure (N m~?) 

Q solution vector (solution variable dependent) 

(r) mean pore radii (m) 

(r?) mean of squared pore radii (m°) 

R discrete residual (kg m~? s7!) 

Rate reaction rate (kmol m~? s7!) 

S source term vector (kg m7? s7!) 

ta anode thickness (m) 

T temperature (K) 

V control volume (m°) 

Xi mole fraction of ith species 

Greek symbols 

A costate variable vector (cost function dependent) 

Qa collision integral 

B design variable vector (design variable depen- 
dent) 

x grid vector (m) 

E Lenard-Jones parameter (m? kg s7?) 

n molecular viscosity (kg ms!) 

v square root of relative molecular weight 

p mass concentration (kg m7?) 

o collision diameter (Å) 

w porosity/totruosity 

Constants 

F Faraday constant 96484.56 (As mol~!) 

K Boltzmann constant 1.3806503 x 107-23 (J K7!) 

Ry universal gas constant 8314.4 (J kmol~! K7!) 

Indices 

d diffusion 

i,j chemical species 

kn Knudsen 

p permeation 

T reforming reaction 


S 


shift reaction 


algorithms to systematically improve existing fuel cell designs, 
particularly when there are many design variables. By using an 
adjoint method, sensitivity derivatives may be obtained for use 
in a design optimization environment. A particular strength of 
adjoint methods is that sensitivity information can be obtained 
with a computational cost that is only weakly dependent on the 
number of design variables, and is therefore enabling technology 
for design studies where many design variables are required. 

In recent years, adjoint methods have been developed and 
utilized for numerical simulations in the aerodynamic com- 
munity for sensitivity analysis, error estimation, and adaptive 
meshing [13—24]. For sensitivity analysis, adjoint methods are 
extremely valuable for determining the derivatives of engi- 
neering cost functions that depend on many design variables. 
Although sensitivity information could theoretically be obtained 
by systematically varying each design parameter independently, 
the cost of repeating the simulation with each variation of a 
design variable renders this approach unusable when a large 
number of parameters are present. The adjoint method is par- 
ticularly suited to this class of problems because the sensitivity 
derivatives can be obtained for all design variables with the com- 
putational cost of a single solution of the nonlinear system used 
for analysis, a single solution of the linear adjoint system, and a 
matrix-vector multiply. 

The primary goal of this study is to formulate and develop 
adjoint methodology for practical applications to fuel cells, with 
particular emphasis on solid-oxide fuel cells. The adjoint method 
is formulated and developed for SOFC applications to provide 
sensitivity derivatives of engineering cost functions. The adjoint 
method is then demonstrated using the one-dimensional model 
of transport phenomena inside an SOFC proposed by Lehnert in 
ref. [5]. Here, sensitivity derivatives with respect to geometric 
and material properties of the anode are obtained for the molar 
concentration of hydrogen at the interface between the anode 
and the electrolyte. Finally, a three-dimensional application is 
included where the cost function is based on the requirement of 
equally distributing fluid through the channels. As a demonstra- 
tion, sensitivity derivatives are obtained for all the mesh points 
defining the surface of the geometry, totaling more than 180,000 
design variables. Grouping some of these design variables to 
describe the shape of the upper and lower walls, the resulting sen- 
sitivity derivatives are then used to reshape the walls to achieve 
improved distribution of fluid amongst the channels. 


2. Discrete adjoint method 


The goal of an adjoint method is to determine sensitivity 
derivatives that can be used in a formal optimization procedure 
for minimizing a specified cost function, which is indicative of 
the performance of the system. In a design setting, a vector of 
design variables is chosen that will be allowed to change during 
the optimization procedure. These design variables may be geo- 
metric in nature or may be indicative of physical parameters such 
as porosity and tortuosity of the anode and cathode or specific 
concentrations of chemical species in the inlet. The sensitivity 
derivatives reflect the changes in the defined cost function with 
respect to each of the design variables. 
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A general optimization procedure begins by first defining a 
meaningful cost function and a desired set of design variables. 
A numerical analysis of the baseline system is then performed. 
The results of the analysis include the solution variables Q of the 
discretized partial differential equations, which are subsequently 
used to determine the initial cost. Because the numerical anal- 
ysis involves discretization of the partial differential equations 
on a computational mesh, it should be noted that Q represents 
the vector of solution variables where each element of the vec- 
tor is representative of one or more physical variables located 
at each mesh point, x. Before reaching a local minimum during 
the optimization procedure, the solution of the partial differential 
equations and the resulting cost function will change in response 
to changes in the design variables. Note that if any of the design 
variables are representative of the shape of the geometry, then 
with each change in these design variables, the mesh must also 
deform to represent the new shape. Ultimately, the cost func- 
tion may have an explicit dependence on the vector of design 
variables, 6, but will also have an implicit dependence because 
Q and x may also depend on the design variables. Therefore, 
the cost function is typically written to indicate the implicit and 
explicit dependence on the design variables as, 


f = f( Q(B), x(B), B) (1) 


After the analysis problem has been solved and the cost 
function has been computed, the next step in the optimization 
procedure is to determine sensitivity derivatives. There are many 
options available for determining the sensitivity derivatives 
including finite differences, forward differentiation [25,26], and 
adjoint methods [13-23]. In this work, adjoint methods are 
used because of their ability to efficiently determine sensitiv- 
ity derivatives for a large number of design variables. Once the 
sensitivity derivatives are obtained, they are used as input to an 
optimization code, such as a quasi-Newton method, to determine 
suitable changes to the design variables to obtain a decrease in 
the cost function. The design variables are updated and the pro- 
cess repeats by analyzing the new configuration that reflects the 
updated design variables. After each analysis, the cost function 
and gradients are typically monitored as an indication of how 
well the design process is proceeding. 


2.1. Sensitivity analysis 


As mentioned in the previous section, the adjoint method is 
used to determine sensitivity derivatives for the cost function, 
f, described in Eq. (1). If R represents the vector of discrete 
residuals at each mesh point, an augmented cost function L can 
be defined in the terms of the original cost function and the vector 
of discrete residuals as following. 


L(Q(B), x(B), B, A= f( Q(B), x(B), B) + ATR(Q(B), x(B), P) 
(2) 


In Eq. (2), A is the vector of Lagrange multipliers (also known 
as costate variables). Note that the augmented cost function, L, 
is a scalar quantity that is identical to the original cost function 
f, when R(Q) is zero, indicating that the steady-state solution is 


obtained. The mesh coordinates and design variables are denoted 
as x and £, respectively. Differentiating the augmented cost 
function with respect to each of the design variables yields the 
following set of equations for dL/08, which is a column vector 
where each element represents the derivative of the augmented 
cost function with respect to a particular design variable. 


aL {af [ax] af a0]! f af [TƏR]! 
Se teak 135+ [ia a} 


dR] fax] far]" 
+ 4/4 A (3) 
op ə | | 9x 
Because the elements of A are arbitrary, the term involving 
the derivatives of the dependent variables with respect to the 
design variables can be eliminated by solving a linear system 


of equations for the costate variables, also known as the adjoint 
equation. 


T 
Palmas) a 
a0 3O 


Once the costate variables are obtained, the derivatives of the 
cost function with respect to all the design variables are obtained 
using a matrix-vector multiplication. 


dL fəf [ax] af aR]? fax] t fart 
ole ATRE + [Bl Lael {4 


(5) 


In numerical simulations, the largest computational cost of 
computing sensitivity derivatives using the adjoint equations is 
due to the solution of the analysis equations and the adjoint 
equation, both of which are independent of the number of design 
variables. The only dependency on the number of design vari- 
ables is in the evaluation of Eq. (5), which is generally much 
cheaper to compute than either the analysis or adjoint solutions. 
The result is a methodology that can be used for computing sen- 
sitivity derivatives for a large number of design variables that is 
only weakly dependent on the number of parameters used in the 
design. 


2.2. Linearization 


Note that the terms in Eqs. (2)-(5) involve differentiation of 
the discrete residual R, the cost function f, and the computa- 
tional mesh x with respect to the dependent variables Q, the 
design variables £, and the location of the mesh points x. Cor- 
rect implementation of this procedure can be extremely tedious 
to accomplish by hand and the resulting code can be difficult 
to maintain. To overcome the difficulties associated with hand 
differentiation, the complex-variable technique of Burdyshaw 
and Anderson [15] and Nielsen and Kleb [23] has been used 
for evaluating all the terms in the matrices required for solving 
the adjoint equations and for evaluating Eq. (5) once the costate 
variables have been obtained. 

The complex-variable technique is derived by expanding 
functions in terms of a complex-variable Taylor series as shown 
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in Eq. (6). This particular example computes the nth element of 
the right hand side vector [d//0Q] of Eq. (4) by adding a complex 
perturbation to the nth element of the Q. 


ih 2 enion 
fo" +im = FQ") + ing (ny + MLO 
(hP f’”"(O") (ih) f(Q”) 

Y 3! t 4! 


h2 n n ht minn 


+ O(h’) 


2! 4! 


h? f(Q”) 


+i nro" BE | +0 © 


where, “i” denotes the imaginary part and “O” is the order of 
accuracy of the next term in the series. The imaginary and real 
parts of the above expression can be written separately as shown 
below. 


h2 Hin ht ENN 
Real [f(Q" + ih)] ~ [ror re i re | 
(7) 
I n ih h2 id n 
ey, [rion - ne J a 


Eqs. (7) and (8) can further be reduced to the following form to 
obtain both the function and the derivative with respect to Q” 
with second-order accuracy. 


f(Q") = Real (f(Q" + ih)) + O(h?) (9) 
f(g”) = == Me +H) 5 oO) (10) 


This technique has been used for evaluating complicated 
algebraic expressions in refs. [27,28] and has been extended for 
cost functions obtained through numerical simulations of par- 
tial differential equations in refs. [11,12]. The significance of the 
equation for the derivative is that by adding a complex pertur- 
bation to the variable of interest and re-evaluating the function, 
sensitivity derivatives can be obtained with relatively modest 
modifications to an existing code. To develop a “forward mode” 
sensitivity analysis code, the primary changes to an analysis code 
involve converting real-valued variables to complex and adding 
a small complex perturbation to the design variable of interest. 
After executing the resulting code, the sensitivity derivatives are 
easily harvested from the complex part of the solution. 

Although the truncation error for this expression has the 
same order of accuracy as a real-valued central finite-difference 
expression, it is important to recognize that the above technique 
does not suffer from accuracy problems due to subtractive can- 
cellation errors, which can often require making adjustments on 
a case-by-case basis when using finite differences. As a result, 
derivatives accurate to numerical precision can be obtained using 
the complex-variable technique by choosing a very small step 
size. The only requirement on the step size is that it should be 
distinguishable from zero as represented by the computer. A sim- 
ple example of using the complex-variable technique is shown 


0.01 T 


Central Finite Difference -$= 
Complex Taylor Series Expansion ---x--- 


sin(x) 
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2 
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Fig. 1. Derivatives obtained using finite-difference and complex-variables. 


in Fig. 1, where the derivatives are obtained for a trigonometric 
function using a real-valued central finite-difference technique 
and the complex-variable technique. Fig. 1 shows the error in the 
computed derivatives obtained using both the central-difference 
and complex-variable methods plotted against the negative of 
the logarithm of the step size. Initially, both methods exhibit 
second-order accuracy as demonstrated by the fact that a one 
order of magnitude reduction in the step size reduces the error 
by two orders of magnitude. As the step size is continually 
reduced, the order of accuracy of the finite-difference method 
eventually degrades and the error begins to increase with sub- 
sequent reduction in the step size. In contrast, the derivatives 
computed using the complex-variable technique continue to 
exhibit true second-order accuracy regardless of the step size. In 
practice, the step size is chosen to be below the square root of 
the machine zero, guaranteeing derivatives accurate to machine 
precision. 

In the context of the present work, the complex-variable tech- 
nique is used to determine the elements of the matrices and 
vectors in Eqs. (2)-(5). In each case, a complex perturbation is 
added to the appropriate variables with respect to which differen- 
tiation is required. For example, to compute the elements of the 
matrix representing the linearization of the discrete residual with 
respect to the flow variables 0R/dQ, a complex perturbation can 
be added to a particular flow variable and the complex-valued 
residual may be computed at each mesh point so that each per- 
turbation effectively yields a column of the matrix. However, 
naive implementation of this technique would be very restric- 
tive for three-dimensional problems due to the large number of 
flow variables. To overcome this difficulty, the method described 
in refs. [15,16] is used to exploit the fact that the computation 
of each residual only depends on mesh points within a finite 
stencil. In this method, the mesh is “colored” into a series of non- 
overlapping stencils to allow multiple columns to be computed 
simultaneously. By following this procedure, the computational 
burden is significantly reduced, requiring only the number of 
residual computations as there are colors in the mesh. Further 
details of this technique are given in refs. [15,16]. A similar 
technique has been described in ref. [23]. 
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A few comments regarding the complex-variable approach in 
comparison to automatic differentiation (AD) techniques such 
as ADIFOR [29], ADIC [30] and TAMC [31] are in order. 
Because the complex-variable method can use arbitrarily small 
step sizes without suffering from subtractive cancellation errors, 
both methods can be used to determine sensitivity derivatives to 
numerical accuracy. In this regard, each method may be viewed 
as an alternative method for achieving the same goals. However, 
each method also has some advantages and disadvantages over 
the other in terms of developing and maintaining software that 
is continually evolving. 

First, each method typically requires some modifications to 
the analysis code before being able to extract derivative infor- 
mation. In the complex-variable approach, the modifications 
are required to change the declaration of floating-point num- 
bers to be of complex data type, to modify comparisons to 
account for mixing of complex and real-valued constants, and 
to sometimes supply intrinsic functions that do not automat- 
ically accept and return complex numbers. In the examples 
given in this paper, no new intrinsic functions had to be sup- 
plied. For automatic differentiation, the analysis code is often 
modified so that the automatic differentiation software can suc- 
cessfully process the analysis source code without errors. In 
this regard, it should be noted that AD software accepts source 
code as input and generates new source code as output that 
includes additional lines of code that provide the capability 
to compute sensitivity derivatives. The types of source-level 
changes required depend somewhat on the computer lan- 
guage and the specific AD software being used. However, the 
changes often include removing non-standard language exten- 
sions, removing arrays of unspecified length, removing common 
blocks from FORTRAN code that are used to pass non-static 
data to subroutines, adding compiler directives to assist in 
the generation of efficient code or preprocessing directives to 
deal with parallelization, eliminating system calls, changing 
function calls to eliminate the use of a single variable for 
both input and output, removing comparisons between incon- 
sistent data types, and removing discontinuous flow control 
[32]. 

A performance comparison of the complex-variable approach 
with automatic differentiation is difficult to state categorically 
because the performance of each approach is dependent on the 
structure of the analysis code, the computer language, compiler 
options and source-level directives passed to the AD software, 
and whether or not further human generated modifications are 
done to the source code generated by the AD software. Experi- 
ence has shown that for the complex-variable approach, software 
written in FORTRAN runs about twice as slow when complex 
arithmetic is used in comparison to when real arithmetic is used. 
For C++, there is a greater variation in run times depending on the 
compiler and the complex classes. The overhead associated with 
using complex arithmetic ranges from about 2—5. Relying on the 
results presented in ref. [32], the range of overhead for automatic 
differentiation is between 2-8, although most of the results are 
within a factor of about 2—4. Based on these results, the complex- 
variable and automatic differentiation approach exhibit roughly 
comparative efficiency. 


One feature of the complex-variable approach is that in the 
process of converting the initial analysis code to accommodate 
complex arithmetic, the code can be made very easy to maintain. 
In particular, as new capability is added or errors are corrected, 
the software can be “self maintaining”. Code generated using 
automatic differentiation software can also be easily maintained 
provided extensive modifications are not required to the AD 
generated source code to obtain efficiency. In this regard, the 
complex-variable approach has the advantage that the source 
code used to obtain sensitivity derivatives is identical to the anal- 
ysis code so it is easy for developers to maintain. Code that has 
been generated from automatic differentiation software includes 
many new lines of code and can be difficult to follow. 

Ultimately, both the complex-variable approach and the use 
of automatic differentiation software are viable solutions for 
obtaining sensitivity derivatives. In this work, the complex- 
variable approach is used because the code is easy to use and 
the code is easy to maintain. 


3. Results 
3.1. One-dimensional SOFC model 


Initial use of the adjoint methodology is demonstrated for 
a numerical simulation similar to that described in ref. [5] 
and illustrated in Fig. 2. Here, the concentration of hydrogen, 
methane, carbon monoxide, carbon dioxide, and water vapor 
are computed through the thickness of an anode with concen- 
trations and pressure prescribed at the boundary of the anode 
with the fuel channel, and the current density is specified at the 
interface with the electrolyte. The simulation assumes a constant 
temperature of 1123 K and includes molecular and Knudsen dif- 
fusion, permeation and reaction kinetics. The anode is assumed 
to be made of standard cermet as described in ref. [5], for which 
properties are given in Table 1. Although details of the govern- 
ing equations can be found in ref. [5], a summary is provided 


Fuel 
Channel 
T 
Anode i ie a 
1 
i 
Electrolyte 
Cathode 
Air 
Channel 


Fig. 2. Computational domain for one-dimensional SOFC anode simulation. 
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Table 1 

Material properties and reaction rate constants from ref. [5] 
kf; 8.0 x 1071! kmol m~? Pa~? s7! 
kb; 1.5 x 10723 kmol m~? Pa~™4 s7! 
kf, 3.2 x 107!! kmol m~? Pa~? s7! 
kbs 3.5 x 107!! kmol m~? Pa~? s7! 
y 1.56 x 107! 

(r) 1.07 x 10-6 m 

0°?) 3.8 x 1071! m 


below to describe the level of physical modeling included in the 
simulations. 

The governing equations used in the simulations are based 
on steady-state mass conservation for each species. 


V- Ji(Q) — Si(Q) = 0 (11) 


In Eq. (11), Q represents the vector of solution variables. For 


; ; T 
this particular problem, Q = [ 0m0, PHz» PCO, PCO); PCH, | y 
where p; denotes the mass concentration of the ith species. 


Ji = MiNi (12) 


e099 
l 


Here, a subscript, “i”, is used to denote each of the species 
included in the simulation and Eq. (12) indicates the relationship 
between mass and molar fluxes. The source term “Sp” indi- 
cates the production and consumption of each species due to 
chemical/electrochemical reactions. 

The molar flux can be subdivided into (1) diffusion flux and 
(2) permeation flux as shown in Eq. (13). 


Ni = NÌ + NP (13) 


The Mean Transport Pore Model (MTPM) [33] is used to com- 
pute the diffusion flux. According to MTPM, 


NS d d 

Nd X ;N; = XiNj 

pet 2 pe = as (14) 
j=l | 
j#i 


where, superscripts “kn” and “m” stand for the Knudsen 
and effective binary diffusion coefficients, respectively. Eqs. 
(15)-(18) illustrate further details involved in computing the 
aforementioned diffusion fluxes [34,35]. 


Dij = wDij (15) 


1 1 
P (an + a) 


Dj; = 1.8823 x 107? (16) 
E poj(2aij 
ae t a (17.1) 
Ce 1.06036 0.193 1.03587 
Vij ~ “F0.1561 " exp(0.476357) | exp(1.529967) 
1.76474 
(17.2) 


+ ——— 
exp(3.894117) 


5 KT 
t (17.3) 
Eij 
Eij = SEF; (17.4) 
2 /8RT\!? 
De = es (eae 1 
i vin (E) (18) 


The permeation flux is given by Eq. (19), and is derived from 
Darcy’s law [33]. 


NP = —X;iB;Vey (19) 


Various terms appearing in this equation can be computed 
using Eqs. (20) and (21). Here, B; is the effective permeability 
coefficient of species i and is described by Eq. (20). Though sev- 
eral different values are suggested for w in the literature [5,33], 
the present simulation uses the value of 7/4. Also, the Knudsen 
number of each component is indicated by Knj;. 


-+ Kn; 2 
B= pet Km.) op (20) 
1+ Kn; 8n 
M;i 
vi = (21.1) 
M mixture 
NS 
Mmixure = X _X;M; (21.2) 


i=l 


The simulation presented in this paper considers two chem- 
ical reactions, namely, methane reforming and water gas shift 
reaction as presented by Eqs. (22.1) and (22.2), respectively. 


kfr 
CH4 + H20€C0O + 3H2 (22.1) 
kb; 
kfs 
CO + H20ÆC02 + H2 (22.2) 


kb, 


Reaction rates are computed using the global reaction model. 
According to this model, 


Rate; = kf, pcn, Po — kbr pco Pix, (23.1) 


Rate, = kf, pco Pmo — kbs pco, Phs (23.2) 


6699 


Subscripts “r” and “s” stand for reforming and shift reactions, 
respectively. Reaction rate constants are taken from [5] and are 
also tabulated in Table 1. The source terms for each species are 
computed using the above reactions rates. 

Two boundary conditions are required to solve Eq. (11) in 
a single spatial direction. Dirichlet boundary conditions are 
applied on the boundary between the anode and the fuel channel, 
i.e. fixed values of mole fractions and pressure are specified. The 
pressure is assumed to be one atmosphere and the mole fractions 
are given in Table 2. At the anode-electrolyte interface the mass 
flux of each species is specified to reflect a prescribed current 
density of 3000 Am~?, which is a Neumann boundary condi- 
tion. The following two electrochemical reactions are accounted 
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Table 2 

Prescribed mole fractions at interface between anode and fuel channel 

Species H? H20 co CO2 CH4 
Mole fraction 0.263 0.493 0.029 0.044 0.171 


for at the anode—electrolyte interface. 


1 
Hp + 502 > H2O (24.1) 


1 
CO + z792 —> CO2 (24.2) 
The current density is related to the species mass flux according 
to Faraday’s law. Thus, conversion rates of the species involved 
in above reactions can be written as the following. 


Ratey, = (—Ratey,0) = (24.3) 


2F 


Rateco = (—Rateco,) = (24.4) 


2F 

To solve the governing equations, the computational domain 
is divided into a finite number of control volumes, and the par- 
tial differential equations are discretized using a finite-volume 
technique. Newton’s method is used to solve the resulting sys- 
tem of nonlinear algebraic equations. Other methods could be 
used with no effect on the final solution. At each iteration, the 
complex-variable technique is used to fill the matrix represent- 
ing the linearization of the discrete residuals with respect to the 
field variables. It should be noted that this matrix is the transpose 
of the matrix required in the development of the adjoint method. 
Although the linearization can be accomplished by hand differ- 
entiation, the complexity of the fluxes shown in Eqs. (13)-(21) 
requires significant effort to obtain an error-free implementa- 
tion. In contrast, linearizing these fluxes is easily accomplished 
using the complex-variable approach. 

For the first simulation, the computational domain is divided 
into 20 equally spaced intervals. Fig. 3 compares the molar for- 
mation rates of all five species from the current simulation with 
those presented by Lehnert et al. [5]. The agreement is satis- 
factory but indicates that there are some modeling differences. 
Specifically, the binary diffusion coefficients and Knudson num- 
bers used in ref. [5] are not given in the reference, and are likely 
different than those employed here. 

To demonstrate the use of the adjoint method for sensitivity 
analysis, the cost function is chosen as the molar concentration 
of hydrogen at the anode-electrolyte interface and the design 
variables include the ratio of porosity to tortuosity, the mean 
pore radii, and the anode thickness. Although only three design 
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Fig. 3. Molar rates of formation through the thickness of the anode. 
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Fig. 4. Distribution of costate variables through thickness of the anode. 


variables are used, many more could be added without signifi- 
cantly impacting the computational cost. Fig. 4 shows the plot 
of costate variables through the anode. Note that the costate 
variables do not typically follow trends exhibited by the species 
formation rates or other field variables such as concentrations, 
but are indicative of the sensitivity of the cost function to changes 
in the residual. This can be seen from Eq. (5), where it is clear 
that in computing the sensitivity derivatives, the costate vari- 
ables weight the contributions due to changes in the residuals 
with respect to the design variables. 


Comparison between sensitivity derivatives obtained using the adjoint method and finite-difference method 


Cost function — molar concentration of hydrogen at anode-electrolyte interface 


ƏLldta aLidw dL/d(r) 
Finite differences 3.0103686566e-01 —1.0500656390e-03 —7.2445001109e + 01 
Adjoint 3.0103704637e-01 —1.0500656734e-03 


—7.244506333e +01 
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Table 3 shows a comparison of sensitivity derivatives 
obtained using the adjoint formulation to those obtained using 
a centered finite-difference technique. For the adjoint results 
shown in Fig. 4, a complex-valued step size of 1 x 1078 is 
used when computing the linearizations. It should be noted that 
smaller step sizes have also been used without changing the 
results, indicating that the derivatives are consistent to machine 
accuracy. The step size for the finite-difference calculations is 
1 x 1075 and has been determined experimentally based on a 
procedure similar to that shown in Fig. 1. As seen, the com- 
parisons of these results are excellent. However, it should be 
noted that accurate sensitivity derivatives were obtained with 
the finite-difference approach by a trial-and-error procedure to 
determine the optimal step size. Secondly, sensitivity deriva- 
tives obtained using finite differences require a solution of the 
governing equations for each design variable. Because this is a 
one-dimensional simulation and there are few design variables, 
this cost is not prohibitive. However, for high-fidelity three- 
dimensional simulations, repeated computations to obtain the 
sensitivity derivatives would be too costly if there are more than 
just a few design variables. 


3.2. Three-dimensional flow problem 


To demonstrate the use of adjoint methodology for a case 
with a large number of design variables, sensitivity derivatives 
are obtained for a three-dimensional fluid flow problem that is 
representative of fuel flowing between channels in a fuel cell. 
The geometry, illustrated in Fig. 5, is similar to that described in 
ref. [36]. For this problem, the fluid enters the cell on the lower 
left and exits on the upper right after flowing through a series of 
manifolds. 

The numerical solution of this problem is obtained by 
assuming incompressible, non-reacting flow and using the three- 


Channel 

mass flow —p> 
measurement 
area 


Fig. 5. Three-dimensional fuel cell geometry. 


Fig. 6. Contours of vertical velocity for initial flow solution. 


dimensional, implicit, unstructured Navier-Stokes flow solver 
described in ref. [37]. This particular geometry is selected as it 
resembles a shape of the planar type SOFC. Note that this sim- 
ulation does not include a high-fidelity physical model required 
for a complete fuel cell simulation. However, the model is use- 
ful for studying general flow distribution between the channels. 
Furthermore, it provides an excellent example for demonstrat- 
ing the use of the adjoint method for determining sensitivity 
derivatives for a large number of design variables. 

The mesh includes approximately 55,000 hexahedral ele- 
ments and about 86,000 nodes. Of these nodes, approximately 
60,000 are distributed along the surfaces of the geometry. Con- 
tours of the computed vertical velocity components are depicted 
in Fig. 6. As seen in the figure, a disparate portion of the mass 
flow traverses through the channels closest to the inlet and exit 
opening. 

To obtain better distribution of fuel amongst the channels, 
a cost function is evaluated along the plane seen in Fig. 5. The 
cost function is the sum of the deviation of the mass flow in each 
channel from the average mass flow through the overall cell [16]. 
In a design optimization, minimizing this cost function would 
improve distribution of the mass flow through each channel. The 
adjoint method has been used to obtain the sensitivity derivatives 
of this cost function with respect to the coordinates of each node 
on the surface of the geometry, yielding a total number of design 
variables in excess of 180,000. Fig. 7 depicts the gradient vectors 
overlain on top of the vertical velocity contours indicating the 
sensitivity of the objective function to changes in the surface 
geometry. Although sensitivity derivatives have been obtained 
for all the surface points in the mesh to demonstrate the utility 
of the adjoint method, an examination of the gradient vectors in 
Fig. 7 suggests that the sensitivity derivatives are largest along 
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Fig. 7. Gradients of design variables. 


the upper and lower walls. New design variables are therefore 
obtained by grouping the original design variables on the upper 
and lower walls to allow ramping of the upper and lower plenum 
surfaces. Similar grouping of design variables in this manner is 
often used to prevent oscillations from appearing in the surface 
and to simplify the design space. Results of the flow distribution 


Fig. 8. Vertical velocity contours after modifying upper and lower walls. 


after one design iteration using steepest descent are shown in 
Fig. 8. As seen in the figure, modifying the shape of the upper and 
lower walls has dramatically improved the distribution of fuel 
through the channels. In addition, the cost function is reduced 
by a factor of four. 


4. Conclusions and future work 


An adjoint based sensitivity analysis technique has been 
developed and demonstrated for applications to solid-oxide fuel 
cells. Sensitivity analysis is demonstrated for a one-dimensional 
transport problem through the anode of a solid-oxide fuel cell, 
where the physical model includes multicomponent diffusion, 
permeation, and appropriate chemical and electrochemical reac- 
tions. The accuracy of sensitivity derivatives is established for 
design variables describing geometric and material properties of 
the anode. 

Because one of the strengths of adjoint methods is in deter- 
mining sensitivity derivatives for large numbers of design 
variables, the adjoint method is further demonstrated on a 
three-dimensional design problem where the objective is to 
achieve equal distribution of fluid between multiple channels 
in a manifold geometry. Sensitivity derivatives are obtained for 
approximately 180,000 design variables describing the geome- 
try of the cell. This problem emphasizes the usefulness of adjoint 
methods in obtaining sensitivity derivatives for a design problem 
with many parameters, which may otherwise be prohibitively 
expensive to compute by individually changing each parameter 
and re-running the entire simulation. The results of the sensitiv- 
ity analysis are used to improve the distribution of fluid amongst 
the channels by grouping the design variables to describe the 
rotational positions of the upper and lower walls. Adjusting the 
walls accordingly results in much improved distribution of fluid 
in the channels. 

Future work is targeted at further developing the adjoint 
method for industrial applications to fuel cell designs. This work 
includes merging more realistic physical models into the three- 
dimensional simulation and adjoint codes as well as coupling 
the sensitivity analysis with formal optimization procedures. 
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